

close all;

%% Load and display the file
struct = load_untouch_nii('img/ET.nii');
cIM = struct.img;


himage = imshow(cIM(:,:,35), []);
% graphical user input for the initial position
p = ginput(1);

% get the pixel position concerning to the current axes coordinates
initPos(1) = round(axes2pix(size(cIM, 2), get(himage, 'XData'), p(2)));
initPos(2) = round(axes2pix(size(cIM, 1), get(himage, 'YData'), p(1)));
initPos(3)=35;

thresVal = double((max(cIM(:)) - min(cIM(:)))) * 0.05;
[nRow, nCol, nSli] = size(cIM);
maxDist = Inf;
tfMean = false;

if initPos(1) < 1 || initPos(2) < 1 ||...
   initPos(1) > nRow || initPos(2) > nCol
    error('Initial position out of bounds, please try again!')
end



simplifyTolerance = 1;

tfSimplify = false;



% initial pixel value
regVal = double(cIM(initPos(1), initPos(2), initPos(3)));

% text output with initial parameters
disp(['RegionGrowing Opening: Initial position (' num2str(initPos(1))...
      '|' num2str(initPos(2)) '|' num2str(initPos(3)) ') with '...
      num2str(regVal) ' as initial pixel value!'])

% preallocate array
J = false(nRow, nCol, nSli);

% add the initial pixel to the queue
queue = [initPos(1), initPos(2), initPos(3)];


%%% START OF REGION GROWING ALGORITHM
while size(queue, 1)
  % the first queue position determines the new values
  xv = queue(1,1);
  yv = queue(1,2);
  zv = queue(1,3);
 
  % .. and delete the first queue position
  queue(1,:) = [];
    
  % check the neighbors for the current position
  for i = -1:1
    for j = -1:1
      for k = -1:1
            
        if xv+i > 0  &&  xv+i <= nRow &&...          % within the x-bounds?
           yv+j > 0  &&  yv+j <= nCol &&...          % within the y-bounds?          
           zv+k > 0  &&  zv+k <= nSli &&...          % within the z-bounds?
           any([i, j, k])       &&...      % i/j/k of (0/0/0) is redundant!
           ~J(xv+i, yv+j, zv+k) &&...          % pixelposition already set?
           sqrt( (xv+i-initPos(1))^2 +...
                 (yv+j-initPos(2))^2 +...
                 (zv+k-initPos(3))^2 ) < maxDist &&...   % within distance?
           cIM(xv+i, yv+j, zv+k) <= (regVal + thresVal) &&...% within range
           cIM(xv+i, yv+j, zv+k) >= (regVal - thresVal) % of the threshold?

           % current pixel is true, if all properties are fullfilled
           J(xv+i, yv+j, zv+k) = true; 

           % add the current pixel to the computation queue (recursive)
           queue(end+1,:) = [xv+i, yv+j, zv+k];

           if tfMean
               regVal = mean(mean(cIM(J > 0))); % --> slow!
           end
        
        end        
      end
    end  
  end
end
%%% END OF REGION GROWING ALGORITHM


struct.img = J;
save_untouch_nii(struct,'img/tt.nii');


